Complexity in parametric Bose-Hubbard Hamiltonians and 
structural analysis of eigenstates 
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We consider a family of chaotic Bose-Hubbard Hamiltonians (BHH) parameterized by the coupling 
strength k between neighboring sites. As k increases the eigenstates undergo changes, reflected in 
the structure of the Local Density of States. We analyze these changes, both numerically and 
analytically, using perturbative and semiclassical methods. Although our focus is on the quantum 
trimer, the presented methodology is applicable for the analysis of longer lattices as well. 
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Understanding the complicated behavior of quantum 
many-body systems of interacting Bosons has been a ma- 
jor challenge for leading research groups over the last 
few years. In fact, the growing theoretical interest was 
further enhanced by recent experimental achievements. 
The most fascinating of these was the realization of Bosc- 
Einstein condensation (BEC) of ultra-cold atoms in pe- 
riodic optical lattices which allows for novel con- 
crete applications of quantum mechanics such as atom 
interferometers and atom lasers. The flood of experi- 
mental realizations includes systems ranging from bond- 
excitations in molecules to cantilever vibrations in 
micro- mechanical arrays jjand Josephson arrays 0. 

The simplest non-trivial model that describes interact- 
ing bosons on a lattice is the Bose-Hubbard Hamiltonian 
(BHH), which incorporates the competition between ki- 
netic and interaction energy of the bosonic system. In 
a substantial part of the existing literature (see for ex- 
ample ||), the dynamics and spectral properties of BHH 
was investigated using a semiclassical picture. In con- 
trast, quantum mechanical calculations of a BHH are of- 
ten limited by severe computational memory restrictions. 
However, it is possible to treat small systems with two or 
three lattice sites. These studies are extremely relevant 
to bond excitations in small molecules few coupled 
Josephson junctions and BECs in optical traps with just 
a few wells Q . As to the experimental realization of the 
last case, microtrap technology || is probably the most 
promising approach for realizing these small systems. Re- 
markably, there is already an experimental realization [i| , 
with promising applications. 

In this context, the two-site system (dimer) has been 
analyzed thoroughly from both the semiclassical [ToL ITU 
and the purely quantum viewpoint [H 0, 0] . Such in- 
vestigations have revealed many interesting phenomena 
like the onset of 7r-phase oscillations, symmetry-breaking, 
and self-trapping of boson population, the latter being 
observed experimentally in Q. In fact, the dimer is in- 
tegrable since there are two conserved quantities, the en- 
ergy and the number of bosons. Nevertheless, the rich- 
ness of the results provides a motivation to go beyond 
the dimer and consider new scenarios where even richer 



dynamics should be observed. The trimer opens new ex- 
citing opportunities, in this respect, since the addition of 
a third site leads to (classically) chaotic behavior, thus 
paving the way to understand longer lattices. The trimer 
has been studied quite extensively in the semiclassical 
regime Jl4. 1151 Il6l|. Surprisingly enough, the quantum 
trimer |l7l Il8j (not to mention longer lattices [l^) is 
barely treated. As a matter of fact, the majority of the 
quantum studies are focused on the statistical properties 
of levels [17]. However, spectral statistics is only the first 
step in understanding the behavior of a complex quantum 
system. Substantially more insight is gained through the 
study of eigenstates. 

In this paper we consider the quantum trimer in the 
chaotic regime and study the structural changes that the 
eigenstates undergo as an experimentally controlled pa- 
rameter k (the coupling strength between neighboring 
sites) is varied. The object of our interest is the over- 
lap of a given perturbed eigenstate \n(ko + 5k)) with the 
eigenstates \m(ko)) of the unperturbed trimer 



P(n\m) = \(n(ko + 6k)\m(k ))\' 



(1) 



Alternatively, if regarded as a function of n for a fixed m, 
the kernel P(n\m) represents (up to some trivial scaling) 
the Local Density of States (LDoS) |2jj. Its lineshape is 
fundamental for the understanding of the associated dy- 
namics since its Fourier transform is the so-called "sur- 
vival probability amplitude" . In our studies, we have 
identified three structural regimes of the P(n\m), which 
are associated with two parametric scales defined as 



Sk qm oc U/N 3/2 and Sk prt cx U/N, 



(2) 



where U = NU, N is the number of interacting bosons 
and U is the on-site boson-boson interaction. For 8k < 
Sk qm the perturbation mixes only neighboring levels: the 
main component of the kernel P(n\m) remains unaffected 
while corrections are captured by standard textbook fi- 
nite order perturbation theory. For Sk qm < 5k < Sk prt a 
non-trivial structure appears, consisting of two distinct 
components: while the tails are still captured by per- 
turbation theory, the central part is of non-perturbativc 
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nature and extends over an energy width T cx N 2 -5k 2 /U . 
For 5k > Sk prt , quantum mechanical perturbation theory 
fails totally. Instead, classical calculations can be used to 
predict the shape of P(n\m). An overview of the para- 
metric evolution of P(n\m) is shown in Fig. 1. 

The trimeric Bosc-Hubbard-Hamiltonian, which de- 
scribes an interacting boson gas confined in a three well 
lattice, is given in second quantization by: 



J ■ 



h = l. 



(3) 



In the BEC framework, k = kg + 5k, is the coupling 
strength between adjacent sites i, j, and can be controlled 
experimentally (in the context of optical lattices this can 
be achieved by adjusting the intensity of the laser beams 
that create the trimeric lattice), while U = A'Kh 2 a s V Bt{ /m 
describes the interaction between two atoms on a single 
site (Kft is the effective mode volume of each site, m is 
the atomic mass, and a s is the s-wave scattering l eng th 
of atoms) . In the context of molecular physics 0, UM k 
represents the electromagnetic and mechanical coupling 
between the bonds of adjacent molecules i,j, while U 
represents the anharmonic softening of the bonds under 
extension. The operators ri, = SjS, count the number 
of bosons at site i; the annihilation and creation opera- 
tors hi and b\ obey the canonical commutation relations 



J i,3 ■ 



Hamiltonian (J3J) has two constants of mo- 



tion, namely the energy E and the total number of par- 
ticles N = J?., n,. Having N = const, implies a finite 
Hilbert-space of dimension M = (JV+2)(JV+l)/2 HIEI- 
An additional 3-fold permutation symmetry allows us to 
reduce further the dimensionality of our space |l7j . 

When N ^> 1 one can adopt a semiclassical point of 
view for Hamiltonian J2J|. Formally, this can be seen 
if we define rescaled creation and annihilation opera- 
tors Cj = 1/yNbi. The corresponding commutators 
[ci,cj] = 5ij/N vanish for N 1 and therefore one can 
treat the rescaled operators as c-numbcrs. The classical 
Hamiltonian H is obtained using the Heisenberg rela- 
tions Ci — > y/Ti exp tVi where ipi is an angle and Ii is the 
associated action. We then get 



NU 



(4) 



The dynamics is obtained from (J3J using the canonical 
equations dli/dt = —dH/dipi and dipi/dt = dH/dli. 
Here t = U ■ t is the rescaled time. The dimensionlcss 
ratio A = k/U @, [H H 03 Gl determines the dy- 
namics of the classical Hamiltonian (@J. For A — ■> the 
interaction term dominates and the system behaves as a 
set of uncoupled sites while for A — > oo the kinetic term 
is the dominant one. In both limits the motion is inte- 
grable. We consider intermediate values of A where the 
classical dynamics is chaotic. The semiclassical limit is 
approached by keeping Aq and U constant while TV — > oo. 



This is crucial in order to keep the underlying classical 
motion unaffected. 

We assume that both Ho = H(Xo) and H = H(X) gen- 
erate classically chaotic dynamics of similar nature. This 
is equivalent with the requirement that 5X = (A— Ao) is 
classically small, i.e., S\ <C 5X c i. An important classical 
observable is the generalized force T(t) = -~(d'H/d\) = 
Y^ijtj \/ lilj exp 1 ^* 3 ^' ~ Vi > . Due to the chaotic dynamics the 
power spectrum C((D) has a finite support fl c \ = 2ir/t erg 
with t eTS being the ergodic time of the system described 
by Eq. (@J. One can use C(ui) as an operative way to 
evaluate 5X C \. The latter is the maximum perturbation 
which keeps C(u>) unaffected. We note that in the exper- 
iment, the relevant parameter that is tuned is the cou- 
pling strength Sk c \ = U ■ 6X c i- We have found that for 

0. 04 < A < 0.2 and in an energy interval H w 0.26 ± 0.02 
the motion is predominantly chaotic. In this regime, and 
for U = 280 we get 5k d = 20 while O c i « 1 (see Fig. 2a). 

Quantum mechanically, we work in the eigenbasis of 
the Hamiltonian Hq. In this basis Hq becomes diagonal 

1. e. E = Em^5 mn where {E$} arc the ordered eigen- 



values. Their mean level spacing A f» 1.517/ N can be 
estimated using TV oc N 2 together with Eq. gj. The 
perturbed Hamiltonian H is written as 



H = E n - 5k B 



(5) 



We mark that although the perturbation strength 5k is 
assumed to be classically small {5k < 5k c \), quantum- 
mechanically it can be very large, i.e., it can mix many 
levels (see the avoided crossings appearing in Fig. 1 for 
large values of 5k). From exact diagonalization we find 
that B is a banded matrix. Its bandprofile can be deter- 
mined using a semiclassical recipe (|B nm | 2 ) « (N 2 /U)A ■ 
C(w nm )/27r [13. The bandwidth is A h = O c i • U. It 
is common to define b = A&/A w 0.6^. The banded 
matrix B and the band profile are illustrated in Fig. 2a. 

A fixed assumption of this work is that 5k <C 5k c \. 
If we require the perturbation to be also quantum 
mechanically small 5k < 5k qm , then we can apply 
standard first order perturbation theory (FOPT) (see 
Fig. 2b). In this case P FOrT (n\m) sa 1 for n = m, 
while P FOPT (n\m) = 5k 2 (|B nm | 2 ) /(E n -E m ) 2 . The per- 
turbation strength <5fc qm is given by the condition that 
only neighboring levels arc mixed, yielding 5k qm = 
A/y/ (|B nm | 2 ). Substituting the expressions for A and 
(|B n m| 2 ) we obtain the results reported in Eq. In 
Fig. 3 we report our numerical results for 5k qm [22]. A 
nice agreement with Eq. J2J is observed. 

For stronger perturbations (5fc qm < 5k one has to em- 
ploy perturbation theory to infinite order. Until now, the 
only formal results regarding Hamiltonians of the type 
(J5J, have been derived by Wigner [2j|. He assumed that 
B nm is a Banded Random Matrix with a flat band profile 
and found that P(n\m) is a Lorentzian 



P plt (n\m) 



5k 2 \B, 



(|B r 



T 2 + {E n ~E m ) 2 ' 



-5k 2 . (6) 
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For T < A Eq. © reduces to P FO p T while for A < T < A b 
the kernel P(n\m) contains two distinct components: a 
central region | E n — E m \ < T where the mixing of levels 
is non-perturbative and a tail region T < \E n — E m \ < At 
which can still be captured by perturbation theory. The 
condition T ~ A& determines <5fc pr t above which the non- 
perturbative core extends all over the bandwidth. There- 
fore © applies as long as Sk < Sk prt oc Sk qm ■ b 5 . In 
the case of non-interacting systems with chaotic classical 
limit recent studies 0, H3 indicated that the above sce- 
nario, based on random matrix modeling, leads to a fairly 
good description of the kernel. Specifically, P(n\m) was 
found to exhibit a core-tail structure where the width 
of the core scaled as T ~ Sk a with 1 < a < 2 |24( in 
contrast to Eq. @. In these studies the parameter T 
was determined (for a given 5k) by imposing normaliza- 
tion of P pr t(n|m). We note that a strict Lorcntzian is an 
idealization of the random matrix modeling. 

Does our BHH model follow the same scenario or will 
the interactions affect the shape of P(n\m)1 Already 
from Fig. 1 we see that as Sk > Sk qm , a non-perturbative 
core starts to evolve and eventually spills over the whole 
bandwidth At. A more detailed comparison (see Fig. 2) 
between P(n\m) and Eq. © shows an excellent agree- 
ment. In the inset of Fig. 3 we report the measured 
width r as a function of Sk. Our numerical data are 
in agreement with a power law behavior T oc 5k a with 
a = 2 ± 0.01. We then investigated the scaling behavior 
of Skp r t |22j . Our numerical data are reported in Fig. 3 
and are in agreement with the results stated in Eq. J5J). 
The latter equation can now be understood in the light of 
Wigncr's arguments since Sk pl -t oc <5fc qm x b 05 and b ~ N. 

For Sk > 6kp T t the core spills over the bandwidth 
and therefore perturbation theory, even to infinite or- 
der, is inapplicable for evaluating P(n\m). In such cases 
one has to rely on completely non-perturbative meth- 
ods 0, H3- For the Wigner model, it was found that 
P(n\m) = l/(27rA)^/4 - ((E n - E m )/A) 2 [23| while in 
systems that have a semiclassical limit the overlap ker- 
nel becomes P(n\m) = tr (p n p m ). Here p TO ({/<}, {<Pi}) 
and p n ({Ii}, {<fi}) are the Wigner functions that cor- 
respond to the eigenstates |m(fco)) and \n(k)) respec- 
tively The trace stands for d{Ii}d{Lpi\ / (2irTi) d inte- 
gration. In the classical limit p can be approximated 
by the corresponding micro-canonical distribution p oc 
S(E— H({Ii}, {ifii})). The latter can be evaluated by pro- 
jecting the dynamics generated by Tto({Ii}, {fi}) = 
onto the Hamiltonian H({Ii} , {<Pi}) = E(t). In Fig. 4a 
we plot the resulting E(t) as a function of time. The clas- 
sical distribution P(n\m) is constructed (Fig. 4b) from 
E(t), averaged over a sufficiently long time. 

In conclusion, we have analyzed the structural changes 
which the eigenstates of a trimcric BHH undergo as the 



coupling strength Sk between the neighboring sites is 
varied. We have found that for Sk < <5fc prt perturba- 
tion theory (to infinite order if Sk > Sk qm ) is applica- 
ble. In this case, the power spectrum of the general- 
ized force C(uj) is an important ingredient for the the- 
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FIG. 1: (Color online) The kernel P(n\m) of the BHH plotted 
as a function of the perturbed energies E n (LDoS represen- 
tation) and for various perturbation strengths Sk > <5fc qm . In 
the lower plane we report the parametric evolution of the en- 
ergy levels within the bandwidth A;, w 230. The energy width 
r is shown (bold line) as a function of 5k. The averaged shape 
of eigenfunctions is given by the same kernel P(n\m) and is 
obtained by just inverting the energy axis. Here, N = 70, and 
Ao = 0.053. 



ory. It is directly experimentally measurable [2(| be- 
cause the momentum distribution of atoms in a lattice 
is ~ y"V exp(ikj)Xj where k is the atomic momentum 

and Xj = ([bj +l bi + h.c]) is the one-particle density ma- 
trix. In the opposite limit Sk > <5fc pr t one can apply 
semiclassical (non-perturbative) considerations. 

It is our pleasure to acknowledge very fruitful discus- 
sions with D. Cohen and G. Kalosakas. 
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